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ABSTRACT 

Four principle models used in cooling water circulation studies 
are introduced. The coupling problems arising when information 
has to be transferred from one model to another are discussed 
and sources of possible errors identified. The errors intro- 
duced when the various equations involved are solved, are 
described together with possible techniques 'to avoid such 
errors. The paper demonstrates that no fail-safe methods are 
available and suggests that results are used only with full 
awareness of the possible errors. 


INTRODUCTION 

In a conference on waste heat management the objectives of the 
modeling work in connection with the design of a large thermal 
plant need no explanation. In this paper we shall look at the 
problems that arise and the possible errors that may occur 
when numerical models are used for predictions of the distri- 
bution of temperatures around the point of outfall and for the 
computation of the transport and dispersion of the waste heat 
in the receiving water body. Estimates of the possibility of 
recirculation and of the ecological impact are based on such 
models. The point that we wish to make is that results of 
such models should be used with a full awareness of problems 
and sources of possible errors. 

Before going into these aspects of numerical modeling it may be 
relevant to set-off numerical modeling against physical model- 
ing. Why IS it that the modeling work is mainly done with the 
use of numerical models? The answer lies in the fact that dis- 
persion of heat and transfer to the atmosphere are difficult 
to model and this is again related to modeling of the different 
scales of turbulence. Abraham - Ref.[l]- demonstrates in fact 
that turbulent stresses and transports can be modeled correctly, 
but only for conditions with a high Reynolds number and where 

Figures referred in the text are given at the end of the paper. 
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a critical value for this number can be established as a solid 
experimental fact. Jet flow is a problem where this criterion 
applies . 

For computations of the transport and dispersion of heat away 
from the outlet - the far-field - mathematical models have to be 
used. The fact that often irregular flow patterns and topogra- 
phies have to be described means that mathematical models usually 
are numerical models rather than analytical. However, here the 
problem of correctly representing the physical processes is 
hardly smaller. This problem is principally related to lack of 
detail in the description in the models. The equations de- 
scribing the instantaneous movements of particles in three 
dimensions are intractable. In order to have a "workable "model 
we have to use time and space averaged forms and this introduces 
dispersion terms. Fickian-type formulations appear to give a 
workable description. However, the wide variety of time and 
length scales involved makes a unified approach with univer- 
sally applicable expressions for the dispersion coefficient in 
these formulations impossible. The correctness of the descrip- 
tion of the physical processes in mathematical form is discussed 
extensively in the literature. Although we may here have a 
first source of errors, we do not intend to discuss these in 
this paper. 


In order to arrive at "workable" models it is furthermore necess 
ary to divide the region in which the processes of dilution, 
transport, dispersion and transfer to the atmosphere take place 
into a near-field around the discharge and a far-field, with 
different models for the processes in these fields. This intro- 
duces the problem of coupling these various models and here 
errors can be introduced. The models that we shall con- 
sider in this connection are a plume model for the near-field 
and hydrodynamic (HD) and transport-dispersion (TD) models for 
the far-field. The plume model we have in mind has the form of 
a Gaussian distribution of velocity and excess temperature around 
a centerline value. The hydrodynamic model in this connection 
is usually an expression of the vertically integrated equations 
of conservation of mass and momentiam over one or two horizontal 
dimensions. Associated with this is a one or two dimensional 
TD-model . For later reference we give here the form for two- 
dimensions 


9c 

9t 


+ 




1 J_ 
h 9x 




where the notation is 


0 ^ 

h h Ax Ay 


( 1 . 1 ) 


c - excess temperature 

CL “ excess temperature discharged 

u,v - horizontal velocity components, integrated over depth. 
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in resp. x- and y-directions 
h - water depth 

Dx/Dy - dispersion coefficients in resp. x- and y-directions 
Ql - discharge of outlet 

a - first order decay factor for heat. 

In order to assess ecological impact the three hydraulic models 
may be supplemented with an ecological mathematical model. 

As will be discussed, coupling problems exist between all four 
models. Apart from errors arising from such problems the 
numerical solution of the actual equations may also introduce 
errors. These are usually related to a necessary discretisation 
of the domain of the models. 

Before continuing with a discussion of possible errors it must 
be remarked that it lies in the nature of the problem that com- 
puted results even when including considerable errors usually 
look "reasonable". The TD-model seldom becomes unstable or 
otherwise indicates its incorrectness . It is difficult to 3 udge 
from a plot of temperature contour lines whether the results are 
correct or incorrect. 


COUPLING PROBLEMS 

Coupling of Near-Field and Far-Field Models 

At the outfall of a thermal power plant a certain mass of water 
with a certain amount of waste heat is discharged into the 
receiving water body with a certain amount of momentum. Using 
a plume model, the distribution of these quantities around the 
outfall can be modeled. The effect of introducing mass and 
momentum should be represented in the HD-model that is used to 
describe the far-field hydrodynamics and the excess temperature 
distribution around the outfall must be transferred to the TD- 
model. However, the resolution in HD and TD-models is usually 
quite different than that used in plume models. In fact most 
HD and TD-models have a discrete representation, whereas plume 
models usually have a continuous representation. In the process 
of transferring quantities computed in the plume model to the 
other models, errors are introduced. 

Coupling Plume and HD-Model 

Buoyancy and remaining jet momentum induce horisontal velocities 
in the receiving water body. If a high velocity surface jet is 
utilized the jet momentum can induce a current pattern which 
can be of considerable importance for the shape of the entire 
temperature field, especially in situations where the currents 
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in the receiving water body are weak. One should consider here 
that the instantaneous shape of the jet usually is a meandering 
plume. Averaged over a certain time the jet appears as a much 
wider 'steady-state' plume. Clearly these effects are difficult 
to represent in, for example, a finite difference HD-model where 
the grid size must be in the order of hundreds of meters for 
reasons of computational economy. 

Coupling HD and TD-model 

The mass of the discharged water can be represented in the HD- 
model as lumped over a few meshes of the grid and this procedure 
needs not introduce errors. However, an error can be introduced 
by incorrect coupling of the equation of mass of the HD-model 
with the equations of the TD-model. The error can be easily 
avoided by following a correct modeling procedure. The point is, 
however, that when this error is made it is not easily detected. 
Another error, that of numerical dispersion, may mask it. 

A HD-model is usually used to create the flow distribution in 
the receiving water body. Neglecting the contribution of the 
source in the mass-equation in this model will hardly be noticed 
in the velocity field obtained as the more important contribution 
of momentum is also - out of necessity - neglected. Observing 
little effect in the HD results may tempt one to neglect the 
source term in the mass equation when deriving the equations for 
the TD-model. We may illustrate this below. 

In deriving the TD-equation (1.1) given in the previous section 
from the principle equation 


9 (he) 
3t 


+ div(hVc) - div(hD grade) + ac = c^ 


( 2 . 1 ) 


, the hydrodynamic mass equation 

at 


+ div(h V ) = ZQl 


( 2 . 2 ) 


, should be used. The source term resulting in eq. 1.1 has the 
form 


Ql(cl~^) 

h Ax Ay 


(2.3) 


However, if the contribution of the source is neglected in eq. 
2.2 the source term will be 


Ql Cl 


h Ax Ay 

The notation here is that of the previous section with further 
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V , the vector (u,v) 

5 , the matrix 

‘-O Dy ** 

The effect on the temperature distribution around the source is 
dramatic; the temperature will in fact rise with each following 
step in time in the computations. The effect is shown in 
Figs, l.a, b and c, where results are presented obtained in 
connection with a site selection study made for a combined de- 
salination and power plant on the east coast of Saudi Arabia. 

The discharge of this plant will be 176 m^/s with an excess 
temperature of 7°C. Fig. l.a shows the most correct solution | 
in which the mass term is included in the equation and where 
no numerical dispersion is present. Fig. l.b shows the same 
situation, now neglecting the mass term, but still without 
numerical dispersion: high temperature around the outlet causes 
the 1°C contour to stretch over a much larger area. Fig. l.c 
shows the same situation now neglecting the mass term and with 
numerical dispersion. One observes that both errors almost 
cancel each other: the results resemble those of Fig. l.a. So 
even with two considerable errors the results may still look 
reasonable! In Fig. l.d results from a hydrodynamic model are 
shown for the strait between Saudi Arabia and the island of 
Bahrain. The velocity field for the TD-model is taken from 
this HD-model. 

Coupling Plume and TD-model 

The mass error and, as we shall see, also numerical dispersion 
can be avoided. Transferring the temperature distribution from 
the plume model to the TD-model presents a more principle 
problem. The usual approach is to assiame that waste heat is dis- 
tributed over a few grid spaces and uniformly mixed from bottom 
to surface. The integral value of excess temperature times 
volume, taken per unit time , is set equal to the waste heat 
discharged, but different volume, excess temperature products 
may be taken to correspond to the same amount of waste heat. 

It IS difficult to make this distribution equivalent to that 
obtained in the plume model where the distribution is principly 
three-dimensional. Using a different niomber of grid spaces 
for the near field distribution gives different far-field con- 
tour lines. The difference is greatest close to the source 
and becomes less further away. However, with a very different 
near-field distribution the picture of the contour lines will 
be very different. An extreme example is given in Figs. 2. a 
and b. The results shown are obtained in connection with the 
site investigation for a 4000 MW nuclear plant in Denmark. 

Fig. 2. a shows the contour picture obtained with a near-field 
distribution over eight grid spacings, whereas in Fig. 2.b 
the distribution is over only one grid spacing. Of course the 
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very high excess temperature in the last case wxll lead to a high 
rate of heat loss to the atmosphere and the contours cover a 
smaller area. Again numerical dispersion could mask the error 
caused by an incorrect near-field distribution. 

The choice of near-field distribution in the TD-model must of 
course be guided by the results obtained in the plume model, 
but because of the different character of the two models the 
transfer will always be imperfect. 

The problem outlined above is, however, not the only principle 
problem in coupling a plume model to a TD-model. The plume 
model assumes a certain temperature of the ambient water. 
Characteristics of the plume such as entrainement and buoyancy 
are dependent on this ambient temperature. One should realise 
in this connection that the 'far-field' is also present near 
the outlet. This is especially the case for a coastal outlet 
with an oscillating tidal current along the coast. In the case 
of a power plant with a large waste heat discharge the tempera- 
ture in the area around the outlet will build up and the 
entraining water will already be heated. This leads us to a 
circular problem in modeling; in order to compute the near-field 
temperature distribution by means of a plume model the ambient 
temperature must be known, but in order to compute the ambient 
temperature one must compute the far-field distribution and this 
is in turn dependent on the near field distribution given by the 
plume. One is forced to spme sort of iterative procedure if 
accurate answers are to be obtained. 

After having looked at these problems of coupling we can con- 
clude that our difficulties stem principally from handling the 
problem in separate models and from differences in resolution 
of these models. This approach is imposed by engineering 
necessity. The problem of difference in resolution may be 
reduced by using a model with varying resolution; detailed around 
the outlet, with a coarser grid away from the outlet. A Finite 
Element Model allows a net with small elements around the outlet 
gradually changing to larger elements further away. In a Finite 
Difference Model a change-of-scale can be used, using a local 
patch of high resolution in an otherwise coarse grid. 

The circular problem of the value of the ambient temperature to 
be used in the pliime computation is hard to overcome. One is 
in fact looking for a "complete" model, where the near-field 
and far-field temperatures are computed simultaneously. Attempts 
have been made to develop such a model using the Marker-in-Cell 
technique. The entire flow region in three dimensions is 
divided into a sufficient number of cells and the computation 
procedure is based on the approximate satisfaction of the inte- 
gral form of the conservation equations for each cell at every 
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time step; the approach xs conceptually that of "box" modeling. 
The demand on computer storage and time is, however, excessive 
and for engineering applications not acceptable. 

Coupling with Ecological Models 

Temperature is an important forcing function in biological and 
chemical processes. Ecological models can be applied to compute 
the consequences of discharges of waste heat. However, such 
models usually deal with slowly varying processes with a time 
scale expressed in weeks and months. The information required 
has for example a form as "number of days per year for which a 
certain temperature level is exceeded" . There is a clear con- 
flict of time scales between the models for temperature distri- 
bution and the ecological models. The approach usually is to 
simulate in the temperature distribution models a short period 
that is statistically equivalent to the period required in the 
ecological model. This last period can be, for example, a year. 
This introduces a statistical interface between the temperature 
distribution models and the ecological model. Errors will be 
introduced, if the quality or quantity of the data on which the 
statistics are to be* based, is insufficient. 

There is also a conflict in the spatial scales between the HD- 
and TD-models and an ecological model. In a complex ecological 
model a large number of ordinary differential equations have to 
be solved for each mesh considered. If a grid of the same mesh 
size as used in the HD- and TD-models would be used in order to 
obtain detailed information on fluxes, the cost of computation 
would become prohibitive. An averaging of hydraulic conditions 
must be introduced with, as a consequence, a loss of accuracy. 

For a more detailed discussion of the use of ecological models 
in power plant site studies we refer to the paper by K.I. Dahl 
Madsen - Ref. [ 2 ]- in this Conference. 


NUMERICAL TECHNIQUES AND ERRORS 

In this section we shall limit the discussion to errors in the 
TD-models and techniques to avoid these. Not that the other 
models have no errors, but particularly in the TD-model the 
errors are difficult to detect and the results may give a false 
impression of correctness. Moreover, the TD-model has a central 
r61e in waste heat studies. 

When the continuous differential equation of transport and dis- 
persion is represented on a discrete grid errors can be intro- 
duced. Well-known is the numerical dispersion that may be intro- 
duced, But although well-known it is difficult to avoid without 
introducing other errors. We may in short recall how this error 
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is introduced. A simple finite difference approximation of the 
terms 



(3.1) 


to 


n+1 



n n n 

c^)/At + u(c^ - Cj_j^)/Ax 


(3.2) 


introduces the truncation error 

u(Ax - uAt) ^ (3.3) 

9x 

, which has the form of a dispersion term. The term depends on 
the choice of the grid and the magnitude of the advective vel- 
ocity. Its value may be many times that of the physical dis- 
persion and may completely invalidate the results . 


In order to avoid this error one may resort to higher order 
difference approximations. However, then a residual numerical 
radiation and skewness appear. The various effects are illus- 
trated in Fig. 3. 


Numerical oscillations also may be introduced by use of an inap- 
propriate scheme. For example, the transport equation 


3c 

3t 


+ 


^c 

3 


0 


(3.4) 


may be approximated by the centred, second order difference 
scheme 



= 0 (3.5) 


However, when initially c = 0 for all j and u is constant, 
provides ^ 
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1 
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A 
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this 


(3.6) 
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which oscillates for all Courant numbers Cr = uAt/Ax, satisfying 

Icr 1<1. 

Higher order approximations have the disadvantage that they 
extend over more and more grid points with increasing orders of 
approximation. This presents problems at the boundaries. Arti- 
ficial assumptions for the approximation have to be made and 
errors are introduced. 

One method to avoid all such types of errors is to apply a 
Lagrangian type of model. The advective term is then in fact 
cut out by moving in a local frame with the local velocity at 
each grid point. When this procedure is followed in a flow 
field varying in time and space, the grid will distort to 
unwieldy forms . Therefore the local frame is moved only for one 
time step and the information obtained is then transferred back 
to points in a fixed grid. One may also express it in a com- 
plimentary form: given a fixed grid at time to + At where was 
the information now in grid point B at time to. The principle 
is illustrated in Fig. 4. Following this principle a practical 
model has been developed at the author's Institute. 

The method requires a sophisticated interpolation technique to 
determine the values of concentration (or temperature) in the 
points at time to# relative to the given fixed grid. The inter- 
polation used is based on a 12-point Everett interpolation. 

The "correctness" of this approach is best illustrated in a so- 
called rotation test in which a Gaussian distribution is rotated 
in a two-dimensional grid around a center point located outside 
the distribution. Results are presented in Fig. 5. It may be 
observed that the shape is fairly well preserved in this rather 
tough test. 

This type of approach may also be used with a Finite Element 
Model. As the F.E.M. technique does not give solutions in dis- 
crete grid points, but as solution surfaces over elements, the 
interpolation is not required (- or is in fact already included 
in the F.E. M. technique) . Clearly higher order elements are 
necessary to obtain results without erroneous dispersion. How- 
ever, compared with finite difference schemes F.E.M. models are 
usually found to be considerably more expensive for time varying 
solutions . 

The 12-point scheme introduced above also requires special for- 
mulations at the boundaries, but satisfactory approximations can 
be obtained. Simply using the correct zero concentration value 
for points beyond a land boundary gives good results. This is 
demonstrated in a so-called L-test in Fig. 6. a. 

The propagation of a circular distribution around a sharp corner 
is depicted for a sequence of time steps. One observes that the 
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circular form is preserved for the stretch before the corner. 
It is dxstorted beyond the corner, but this is physically re- 
alistic. The deformation is caused by the flow field around 
the corner. In order to illustrate how an incorrect scheme 
with considerable numerical dispersion may distort a circular 
distribution, corresponding results of such a scheme are shown 
in Fig. 6.b. 


At open water boundaries also special attention is required. 
However, satisfactory results are usually obtained when such 
boundaries are recognized as being either inflow or outflow 
boundaries, with an assiimption on the mixing conditions in the 
water body beyond the boundary. 


A quite different method, that probably is the most accurate, 
is based on a spectral technique. The method is developed for 
an application in atmospheric pollution by Christensen and 
Pram - Ref . 4 . The technique was applied to hydraulic problems 
in the author's institute. The method is very accurate, but 
computer costs are about four t'imes as high as for the 12-point 
scheme. The method also has limitations with regard to resol- 
ution of realistic topographies . In short the technique is as 
follows: It is assumed that c can be approximated by a finite 

Fourier representation 

^ ikx 

f(x,t) = I A(k,t) e^ , with i = (3.8) 

k 


For a given continuous function f(x,t) one can always find 
A(k,t) such that f(x,t) = c(x,t) on grid points x^ = jAx: 

A(k,t) = i ? c (t) 

The simple equation 


(3.9) 


3c 

3t 


+ 



0 


. 10 ) 


transforms to 


3A(k,t) 

3t 


u i k A (k , t) = 0 


(3.11) 


From this equation A(k,t) can be computed at each new time nAt 
with A(k,t) given, c can be obtained from an inverse Fourier- 
transform. The point here is to note that the advection term 
in (3.11) IS not approximated by a finite difference form so 
that numerical dispersion is avoided. 
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DISCUSSION 

A discussion of modeling problems and errors may leave the reader 
with a rather gloomy impression of the whole modeling effort. 

And there is more, in addition to the problems discussed there 
are the difficulties encountered when dispersion coefficient and 
heat transfer coefficients have to be selected. Field investi- 
gations before the plant has been constructed can only give a 
very limited impression of the mixing characteristics of the 
receiving water body, as discharge volume, discharge momentum 
and buoyancy cannot be represented. 

However, if the models are used with an awareness of the inac- 
curacies and if sensitivity tests are made for the important 
parameters, usefull predictions can be made. Such predictions 
would allow design considerations of outlet-intake construction 
cost against possible recirculation as discussed in the paper 
by Schr0der - Ref. [5]- in this Conference, and an evaluation of 
ecological consequences . 

We may underline the above statement with a final result. The 
possibility of simulating the transport-dispersion process over 
an extended period of time in a realistic topography is illus- 
trated in Fig. 7. The transport and dispersion of a conservative 
substance with an irregular release - 16 hrs out of 24 per day - 
IS simulated. The area concerned is K0ge Bay, south of Copen- 
hagen. In a sequence of plots the development of the concen- 
tration field from an initial distribution to the situation 
after one week is shown. The results after one week are com- 
pared to measurements. One may observe that after having been 
through considerable changes, the characteristics of the field 
after seven days compare well with measurements. 
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Fig. 3: Examples of Errors in Fig. 4: Principle of 12- 

Difference Approxi- Point Scheme, 

mations to Transport- 
Dispersion Equations. 
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(e) - computed 168 hrs (f) - recorded 168 hrs 

Fig. 7. e-f: Sequence of Recorded and Computed Concentration 

Contours for K0ge Bay, Continued. 



(g) - typical flow pattern at (h) - typical flow pattern at 

27 hrs, southgoing 51 hrs, northgoing 

current in the Sound, current in the Sound, 

westerly winds 10-15 m/s southerly winds 8-10 m/s 


Fig. 7. g-h; Typical Flow Patterns. 




